source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/SLE_TRPM2_MfMo/'
kn.pal2 <- c('#d44','#485ffc')

# GSE125527 ---------
platform_meta <- GEOquery::getGEO('GSE125527')

platform_meta$`GSE125527-GPL20301_series_matrix.txt.gz` |>
  pData() |>
  as_tibble() |>
  pull(title)

platform_meta$`GSE125527-GPL24676_series_matrix.txt.gz` |>
  pData() |>
  as_tibble() |>
  pull(title)

GEOquery::getGEOSuppFiles('GSE125527', fetch_files = F)

ptid <- read_csv('https://ftp.ncbi.nlm.nih.gov/geo/series/GSE125nnn/GSE125527/suppl/GSE125527%5FoldPatientId%2DnewPatientId.csv.gz')

ptid

idcol <- read_csv('https://ftp.ncbi.nlm.nih.gov/geo/series/GSE125nnn/GSE125527/suppl/GSE125527%5Fcell%5Fid%5Fcolnames.csv.gz')

idcol

gene_csv <- read_csv('https://ftp.ncbi.nlm.nih.gov/geo/series/GSE125nnn/GSE125527/suppl/GSE125527%5Fgene%5Fid%5Frownames.csv.gz', col_names = 'gene')

gene_csv |>
  filter(gene %in% c('TRPM2','KNP3','LTRPC2','NUDT9L1','NUDT9H','EREG1'))

gene_csv |>
  filter(str_detect(gene, 'TRP'))

cellmeta <-
  read_csv('https://ftp.ncbi.nlm.nih.gov/geo/series/GSE125nnn/GSE125527/suppl/GSE125527%5Fcell%5Fmetadata.csv.gz',
           col_names = c('id','barcode','sample','tissue','group','celltype'),
           skip = 1)

cellmeta |> dplyr::count(celltype)

## raw data ---------
boland_ena <- read_tsv('~/append-ssd/alaria2/boland20uc/PRJNA516681.tsv')

boland_ena |>
  filter(str_detect(sample_title, 'PBMC_10x_scRNA')) |>
  write_tsv('~/append-ssd/alaria2/boland20uc/boland20rna_pbmc.tsv')

boland_ena |>
  filter(str_detect(sample_title, 'R_10x_scRNA')) |>
  write_tsv('~/append-ssd/alaria2/boland20uc/boland20rna_rectum.tsv')

### nxf input ---------
nxf_pbmc <- read_csv('~/append-ssd/alaria2/boland20uc/pbmc/nxf_input.csv')

nxf_pbmc |>
  mutate(sample = str_c(sample, '_', seq_along(sample)), .by = readgroup) |>
  write_csv('~/append-ssd/alaria2/boland20uc/pbmc/nxf_input.csv')

nxf_rec <- read_csv('~/append-ssd/alaria2/boland20uc/rectum/nxf_input.csv')

### nxf output --------
nxfo_pbmc <- list.dirs('/home/supervisor/mist2/gjsx/boland20pbmc', recursive = T) |>
  str_subset('filtered_full$')

pbmc_mex <-
nxfo_pbmc |>
  str_extract('(C|U).+scRNA_\\d') |>
  str_remove('10x_scRNA_') |>
  str_replace_all('_','.') |>
  set_names(x = nxfo_pbmc, nm = _) |>
  Read10X()

pbmc_mex |> glimpse()

sobj_pbmc |>
  GetAssayData() |>
  DropletUtils::write10xCounts(path = 'boland20pbmc_scRNA.h5', x = _)

pbmc_mex <- Read10X_h5('boland20pbmc_scRNA.h5')

foo <- pbmc_mex[1:5,1:5]

foo

n_cell <- tibble(barcode = colnames(pbmc_mex)) |>
  count(barcode)

n_cell |>
  filter(n > 1) |>
  mutate(sample = str_remove(barcode, '.{17}$')) |>
  ggplot(aes(sample)) + geom_bar()

n_cell |>
  filter(n == 1) |>
  mutate(sample = str_remove(barcode, '.{17}$')) |>
  ggplot(aes(y = sample)) + geom_bar()

dedup_cell <- n_cell |>
  filter(n == 1) |>
  pull(barcode)

pbmc_mex <- pbmc_mex[,dedup_cell]

### rectum star ---------
nxfo_rect <- list.dirs('/home/supervisor/mist2/gjsx/boland20rectum/', recursive = T) |>
  str_subset('filtered_full$')

rect_mex <- nxfo_rect |>
  str_extract('(C|U).+scRNA_\\d') |>
  str_remove('10x_scRNA_') |>
  str_replace_all('_','.') |>
  set_names(x = nxfo_rect, nm = _) |>
  Read10X()

rect_mex |> glimpse()

rect_mex <- Read10X_h5('boland20rectum_scRNA.h5')

n_cell <- tibble(barcode = colnames(rect_mex)) |>
  count(barcode)

n_cell |>
  filter(n > 1) |>
  mutate(sample = str_remove(barcode, '.{17}$')) |>
  ggplot(aes(sample)) + geom_bar()

n_cell |>
  filter(n == 1) |>
  mutate(sample = str_remove(barcode, '.{17}$')) |>
  ggplot(aes(y = sample)) + geom_bar()

dedup_cell <- n_cell |>
  filter(n == 1) |>
  pull(barcode)

rect_mex <- rect_mex[,dedup_cell]

## PBMC ---------
sobj_pbmc <- pbmc_mex |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj_pbmc <- sobj_pbmc |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj_pbmc |> VlnPlot('mito_ratio', pt.size = 0)

sobj_pbmc |> dim()

# 507k
sobj_pbmc <- sobj_pbmc |>
  filter(mito_ratio < 10, !(orig.ident %in% c('C17.PBMC.1', 'C18.PBMC.1')))

# 197k
sobj <- sobj_pbmc |>
  filter(nFeature_RNA > 500)

sobj_pbmc |>
  ggplot(aes(nFeature_RNA)) + geom_density() +
  geom_vline(xintercept = c(400,500))

sobj <- sobj |>
  mutate(sample = str_extract(orig.ident, '(C|U)\\d+'),
         group = ifelse(str_detect(sample, '^C'), 'HC', 'UC'))

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident','sample'))

sobj |> DotPlot(seurat_markers, cols = 'RdYlBu', cluster.idents = T)

sobj <- sobj |>
  mutate(celltype = case_match(as.numeric(seurat_clusters),
                               c(14,6,19)~'B cell',
                               20 ~ 'Plasma cell',
                               c(18,12,4)~'CD8 T cell',
                               5~'NK cell',
                               c(2,9,3,17,8)~'CD4 T cell',
                               c(21,13,15) ~ 'DC',
                               .default = 'Monocyte'
  ))

myeloid_markers <- list(
  Macrophage = c("CD68", "CD163", "MSR1"),
  Monocyte = c("CD14", "FCGR3A", "S100A8", "C1QA"),
  Dendritic = c("CD1C", "CLEC9A", "CD83"),
  Neutrophil = c("CSF3R", "FCGR3B", "S100A9"),
  pDC = c('LILRA4','CLEC4C','TCF4'),
  'TRPM2'
)

sobj |>
  DotPlot(myeloid_markers, cols = 'RdYlBu', cluster.idents = T) +
  RotatedAxis()

sobj |>
  DimPlot(group.by = 'celltype', cols = 'Paired')

### TRPM2 dotplot in all ---------
m2_2d <- sobj |>
  DotPlot2d('TRPM2', group.x = group, group.y = celltype) |>
  pluck('data')

m2_2d |>
  mutate(group.y = fct_reorder(group.y, avg.exp)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type',
       title = 'TRPM2 expression in UC patient PBMC') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('UC.PBMC.TRPM2.dotplot')

### M2 featureplot ---------
sobj |> FeaturePlot('TRPM2', cols = c('lightgrey','red'), order = T,
                    split.by = 'group', raster = F) &
  theme_jpub

### rds save ------------
sobj |> write_rds('mission/SLE_TRPM2_MfMo/boland20pbmc.rds')

sobj_mo <-
  sobj |> filter(celltype == 'Monocyte')

sobj_mo |> write_rds('mission/SLE_TRPM2_MfMo/boland20pbmc_Mo.rds')

### M2 expr by sample ---------
sobj_mo |> DotPlot('TRPM2', cols = 'RdYlBu', group.by = 'orig.ident')

m2_by_sample <- last_plot() |>
  pluck('data')

m2_by_sample |>
  mutate(group = ifelse(str_detect(id, '^U'), 'UC', 'HC')) |>
  ggplot(aes(group, avg.exp, color = group)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .1, height = 0) +
  theme_pubr() +
  stat_compare_means(method = 't.test', comparisons = list(c('HC','UC')),
                     color = 'black') +
  labs(x = 'Group', y = 'Average expression',
       title = 'TRPM2 expression in ulcerative colitis PBMC')
  
### MF M2-dotplot & UMAP ---------
sobj_mo <- sobj_mo |> quick_process_seurat(skip_norm = T, res = .4)

sobj_mo |> DimPlot(cols = paired.pal, label = T, label.box = T, repel = T)

sobj_mo |>
  DotPlot2d('TRPM2', seurat_clusters, group)

sobj_mf |>
  filter(group == 'IBD') |>
  DotPlot('TRPM2', cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Macrophage clusters')

mf_leiden_m2 <- last_plot() |>
  pluck('data')

mf_leiden_m2 |>
  mutate(seurat_clusters = id, avg.exp.scaled, avg.exp, .keep = 'none') |>
  left_join(x = sobj_mf, y = _) |>
  ggplot(aes(umap_1, umap_2, color = avg.exp.scaled)) +
  geom_point(size = 1) +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_classic() +
  labs(color = 'Mean expr',
       title = 'Mean expression of TRPM2 in macrophages clusters')

## rectum -----------
sobj_rect <- rect_mex |>
  CreateSeuratObject(min.cells = 3, min.features = 400)

sobj <- sobj_rect |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj |> dim()

# 191k
sobj <- sobj |>
  filter(mito_ratio < 10, !(orig.ident %in% c('C17.R.1', 'C18.R.1')))

sobj |>
  ggplot(aes(nFeature_RNA)) + geom_density() +
  geom_vline(xintercept = c(500))

sobj <- sobj |>
  mutate(sample = str_extract(orig.ident, '(C|U)\\d+'),
         group = ifelse(str_detect(sample, '^C'), 'HC', 'UC'))

sobj <- sobj |>
  quick_process_seurat(batch = c('orig.ident','sample'))

sobj |> DotPlot(seurat_markers, cols = 'RdYlBu', cluster.idents = T)

sobj <- sobj |>
  mutate(celltype = case_match(as.numeric(seurat_clusters),
                               c(14,6,19)~'B cell',
                               20 ~ 'Plasma cell',
                               c(18,12,4)~'CD8 T cell',
                               5~'NK cell',
                               c(2,9,3,17,8)~'CD4 T cell',
                               c(21,13,15) ~ 'DC',
                               .default = 'Monocyte'
  ))

myeloid_markers <- list(
  Macrophage = c("CD68", "CD163", "MSR1"),
  Monocyte = c("CD14", "FCGR3A", "S100A8", "C1QA"),
  Dendritic = c("CD1C", "CLEC9A", "CD83"),
  Neutrophil = c("CSF3R", "FCGR3B", "S100A9"),
  pDC = c('LILRA4','CLEC4C','TCF4'),
  'TRPM2'
)

sobj |>
  DotPlot(myeloid_markers, cols = 'RdYlBu', cluster.idents = T) +
  RotatedAxis()

sobj |>
  DimPlot(group.by = 'celltype', cols = 'Paired')

### TRPM2 dotplot in all ---------
m2_2d <- sobj |>
  DotPlot2d('TRPM2', group.x = group, group.y = seurat_clusters) |>
  pluck('data')

m2_2d |>
  mutate(group.y = fct_reorder(group.y, avg.exp)) |>
  BubblePlot(d2 = T) +
  labs(x = 'Group', y = 'Cell type',
       title = 'TRPM2 expression in UC patient PBMC') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_source_plot('UC.PBMC.TRPM2.dotplot')

### M2 featureplot ---------
sobj |> FeaturePlot('TRPM2', cols = c('lightgrey','red'), order = T,
                    split.by = 'group', raster = F) &
  theme_jpub

### rds save ------------
sobj |> write_rds('mission/SLE_TRPM2_MfMo/boland20rectum.rds')

sobj_mf <-
  sobj |> filter(seurat_clusters == 17)

sobj_mf |> write_rds('mission/SLE_TRPM2_MfMo/boland20rectum_MF.rds')

### M2 expr by sample ---------
sobj_mf |> DotPlot('TRPM2', cols = 'RdYlBu', group.by = 'orig.ident')

m2_by_sample <- last_plot() |>
  pluck('data')

m2_by_sample |>
  mutate(group = ifelse(str_detect(id, '^U'), 'UC', 'HC')) |>
  ggplot(aes(group, avg.exp, color = group)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .1, height = 0) +
  theme_pubr() +
  stat_compare_means(method = 't.test', comparisons = list(c('HC','UC')),
                     color = 'black') +
  labs(x = 'Group', y = 'Average expression',
       title = 'TRPM2 expression in ulcerative colitis rectum macrophage')

### MF M2-dotplot & UMAP ---------
sobj_mo <- sobj_mo |> quick_process_seurat(skip_norm = T, res = .4)

sobj_mo |> DimPlot(cols = paired.pal, label = T, label.box = T, repel = T)

sobj_mo |>
  DotPlot2d('TRPM2', seurat_clusters, group)

sobj_mf |>
  filter(group == 'IBD') |>
  DotPlot('TRPM2', cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Macrophage clusters')

mf_leiden_m2 <- last_plot() |>
  pluck('data')

mf_leiden_m2 |>
  mutate(seurat_clusters = id, avg.exp.scaled, avg.exp, .keep = 'none') |>
  left_join(x = sobj_mf, y = _) |>
  ggplot(aes(umap_1, umap_2, color = avg.exp.scaled)) +
  geom_point(size = 1) +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_classic() +
  labs(color = 'Mean expr',
       title = 'Mean expression of TRPM2 in macrophages clusters')